/*
 * morphologyOperations.cpp
 *
 *  Created on: Dec 15, 2013
 *      Author: justin
 */

#include "morphologyOperations.h"
#include "pbStatements.h"
#include <iostream>
using namespace std;
int child_row_cmp(void* _a,void* _b,void* avgheight)
{
    childImg* a =*(childImg**)_a;
    childImg* b =*(childImg**)_b;
    return (abs(a->y - b->y) < *(float *)avgheight) ? (a->x - b->x) : a->y - b->y;
}
int coor_row_cmp(void* _a,void* _b,void* avgheight)
{
    CvPoint* a =(CvPoint *)_a;
    CvPoint* b =(CvPoint *)_b;
    return (abs(a->y - b->y) < *(float *)avgheight) ? (a->x - b->x) : a->y - b->y;
}
void rotateContour(CvSeq * contour, float angle, int centerX, int centerY, bool isInner) {
    if(angle == 0) return;
    int tot = contour->total;
    CvPoint pt;
    float vcos = cos(angle);
    float vsin = sin(angle);
    for(int i=0;i<tot;i++) {
        cvSeqPop(contour,&pt);
        //旋转公式：
        //x'= (x - rx0)*cos(RotaryAngle) + (y - ry0)*sin(RotaryAngle) + rx0
        //y'= =-(x - rx0)*sin(RotaryAngle) + (y - ry0)*cos(RotaryAngle) + ry0
        pt.x = (float)(pt.x - centerX)*vcos + (float)(pt.y - centerY)*vsin + centerX;
        pt.y = -(float)(pt.x - centerX)*vsin + (float)(pt.y - centerY)*vcos + centerY;
        cvSeqPushFront(contour,&pt);
    }
    if(contour->v_next) {
        rotateContour(contour->v_next, angle, centerX, centerY,1);
    }
    if(isInner && contour->h_next) {
        rotateContour(contour->h_next, angle, centerX, centerY,1);
    }
}
float inclineDetect(CvSeq * coors, int avgHeight) {
    int periorY = -1;
    int tot = 0;
    float k = 0, r = 0;
    float xya = 0, xa = 0, ya = 0, xsa = 0, ysa = 0;
    for(int i=coors->total;i>0;i--) {
        CvPoint pt;
        cvSeqPopFront(coors,&pt);
        if(periorY != -1 && abs(pt.y - periorY) > avgHeight) {
            xya /= tot;
            xa /= tot;
            ya /= tot;
            xsa /= tot;
            ysa /= tot;
            float kt = (xya-xa*ya)/(xsa-xa*xa);
            float rt = abs((xya-xa*ya)/sqrt((xa*xa-xsa)*(ya*ya-ysa)));
            cout<<kt<<", "<<rt<<endl;
            if(rt > r) k = kt;
            r = rt;
            tot = 0, xa = 0, ya = 0, xsa = 0, ysa = 0, xya = 0;
        }
        xya += pt.x*pt.y;
        xa += pt.x;
        ya +=pt.y;
        xsa += pt.x*pt.x;
        ysa += pt.y*pt.y;
        periorY = pt.y;
        tot++;
    }
    if(abs(atan(k)) > 0.01) return atan(k);
    else return 0;
}
void morphologyOperations(IplImage * rawImage, CvSeq * &childs, CvMemStorage * &stofch) {
    IplImage *dst = cvCreateImage(cvGetSize(rawImage), 8, 1);
    cvCvtColor(rawImage, dst, CV_RGB2GRAY);
    cvSmooth(dst, dst, CV_GAUSSIAN, 1);
    cvThreshold(dst, dst, 0 , 255, CV_THRESH_OTSU);
    cvShowImage("a", dst);
    CvMemStorage *storage=cvCreateMemStorage(0);
    CvMemStorage *stofcoor=cvCreateMemStorage(0);
    stofch = cvCreateMemStorage(0);
    CvSeq * contours=0;
    cvFindContours(dst,storage,&contours,sizeof(CvContour), CV_RETR_CCOMP, CV_CHAIN_APPROX_NONE);
    CvSeq *p = contours;
    childs = cvCreateSeq(0, sizeof(CvSeq), sizeof(childImg *), stofch);
    CvSeq *coors = cvCreateSeq(0, sizeof(CvSeq), sizeof(CvPoint), stofcoor);
    float avgHeight = 0;
    for (;p!=0;p=p->h_next) {
        double tmparea = fabs(cvContourArea(p));
        for(CvSeq * pt = p->v_next;pt;pt=pt->h_next) tmparea -= fabs(cvContourArea(pt));
        CvRect aRect;
        aRect = cvBoundingRect(p, 1);
        float whratio = (float)aRect.width/aRect.height;
        int boxsize = aRect.width * aRect.height;
        if(tmparea < 20 || whratio > 1 || whratio < 0.5 || tmparea/boxsize > 0.7) {
            CvSeq *pOldSeq = p;
            if(p->v_next) {
                cvClearSeq(p->v_next);
            }
            if(pOldSeq->h_prev)
            {
                p = pOldSeq->h_prev;
                p->h_next = pOldSeq->h_next;
                if(pOldSeq->h_next) pOldSeq->h_next->h_prev = p;
                cvClearSeq( pOldSeq );
            }
            else
            {
                p->h_next->h_prev = NULL;
                contours = p->h_next;
                cvClearSeq( pOldSeq );
            }
            continue;
        }
        avgHeight += aRect.height;
        CvPoint pt = cvPoint(aRect.x+aRect.width/2, aRect.y+aRect.height/2);
        cvSeqPushFront(coors,&pt);
    }
    avgHeight /= coors->total;
    cvSeqSort(coors, (CvCmpFunc)coor_row_cmp, &avgHeight);
    float incline = inclineDetect(coors, avgHeight);
    for(p=contours;p!=0;p=p->h_next) {
        rotateContour(p, incline, cvGetSize(rawImage).width/2, cvGetSize(rawImage).height/2,0);
    }
    for(p=contours;p!=0;p=p->h_next) {
        CvRect * aRect = new CvRect;
        *aRect = cvBoundingRect(p, 1);
        IplImage * chimg = cvCreateImage(cvSize(aRect->width, aRect->height), 8, 1);
        childImg * child = new childImg(chimg, aRect->x, aRect->y, aRect->width, aRect->height);
        for (int i = 0; i < cvGetSize(chimg).height; i++)
        {
            for(int j = 0; j < cvGetSize(chimg).width; j++)
            {
                ((uchar *)(chimg->imageData +  i * chimg->widthStep ))[j] = 255;  //对指定像素进行修改
            }
        }
        cvSeqPushFront(childs, &child);
        cvDrawContours(chimg, p, cvScalar(0), cvScalar(0), -1, CV_FILLED, 8, cvPoint(-aRect->x, -aRect->y));
    }
    cvSeqSort(childs,(CvCmpFunc)child_row_cmp, &avgHeight);
    cvReleaseMemStorage(&storage);
    cvReleaseImage(&dst);
}

